Journal of Computational Chemistry

Protein recognition mimicry is of great interest in the field of molecular bioengineering and rational design, with mutagenesis frequently employed to analyze the effects of altering amino acids involved in molecular recognition. The conformational and energetic effects of such alterations can be investigated in detail with the help of molecular dynamics (MD) methodologies. While existing MD-based computational tools can be used to explore a particular set of mutations at a time, suitable for small-scale studies, high-throughput (HTP) exploration of protein recognition for engineering purposes would greatly benefit from an integrative platform that streamlines preparation, mutagenesis, simulation and post-processing of up to several thousand molecular systems, along with robust tools for comprehensive and straightforward comparative analysis. DyME (Dynamic Mutagenesis Engine) is a distributed platform that enables systematic investigations of protein recognition mimicry by combining HTP mutagenesis, solvated MD simulations and a Toolbox for comparative analysis (TCA), including interfacial water-site mapping. DyME uses 3D structural information of any protein-protein or protein-DNA complex as input. Its automated MD-based mutagenesis engine facilitates systematic investigation of how site-specific alterations affect recognition, enabling the organization of single, double and triple modifications into combinatorial libraries for comprehensive comparative analysis. In DyME, relevant MD trajectory-derived data is scavenged and stored into a central database, providing aggregation capabilities that ease multi-feature analysis across an extensive collection of simulations. An interactive web-GUI and specialized widgets simplify preparation and efficient molecular and numerical comparative exploration. DyMEs capabilities are evaluated using available experimental data. Its source code is available at https://github.com/pisabarro-group/DYME

Recently, a number of tools have been released that generate ensembles of protein structures based on artificial intelligence (AI) approaches. Although ensembles generated by the tools differ significantly, we demonstrate a computational path to harmonizing the various outputs under a stationary condition using two complementary physics-based approaches. In the first stage, the AI ensemble is used to seed a weighted ensemble (WE) simulation, promoting relaxation toward the steady state. In the second stage, trajectory segments generated by WE are reweighted to steady state using the recently developed RiteWeight (RW) algorithm. We applied this approach to generate an atomically-detailed equilibrium ensemble of unliganded adenylate kinase conformations, starting from ensembles produced by three AI tools: AFSample2, ESMFlow-PDB (trained from PDB structures), and ESMFlow-MD (trained from molecular dynamics simulation data). Dramatic differences in the AI-generated ensembles are largely erased during the WE-RW process, yielding a consistent description of the equilibrium ensemble for a given force field.

The incorporation of non-canonical amino acids (ncAAs) is a powerful strategy for introducing novel chemical functions into proteins. Molecular dynamics (MD) simulations are essential for understanding the structural and dynamic effects of these modifications, yet the creation of accurate force field parameters for ncAAs remains a significant bottleneck. Current parameterisation methods are often inaccurate or computationally expensive. To address this, we present drFrankenstein, an automated pipeline for generating AMBER force field parameters for ncAAs. drFrankenstein is a robust and accessible tool that streamlines the parameterisation workflow, enabling the routine use of MD simulations to study the behaviour of ncAA-containing proteins.

There is a growing demand for molecular dynamics simulations to explore longer timescale behavior of giant protein-DNA complexes such as chromatin. To address this need, we extended OpenCafeMol, a GPU-accelerated residue-level coarse-grained molecular dynamics simulator originally developed for proteins and lipids, to support 3SPN.2 and 3SPN.2C DNA models. We also implemented a hydrogen-bond-type many-body potential to model DNA-protein interactions more accurately. To further improve computational efficiency, we introduced a localized scheme for calculating base-pairing and cross-stacking interactions. Benchmark tests show that OpenCafeMol on a single GPU achieves up to 200-fold speed-up for DNA-only systems and up to 100-fold speed-up for DNA-protein complexes compared to CPU-based simulations. To demonstrate the capability of our implementation for long-timescale biological processes, we simulated an archaeal SMC-ScpA complex undergoing DNA translocation via segment capture (a proposed mechanism for DNA loop extrusion) in the presence of a DNA-bound obstacle. We observed continuous captured-loop growth accompanied by obstacle bypass within the segment capture framework.

Wild-type T4 lysozyme (T4L) is used as a benchmark to evaluate conformational sampling across generative AI, AI-accelerated molecular simulation (AMS), and physics-based enhanced molecular dynamics (EMD). A four-state model: exposed/open, exposed/closed, buried/open, and buried/closed; is defined using physically meaningful collective variables. While generative AI methods (AF-cluster, MSA subsampling of AlphaFold2, ConforFold, AlphaFlow, ESMFlow, ConfRover, BioEmu) largely sample only the exposed/open state, AMS integrating generative ensembles with iterative molecular dynamics, recovering all states and reproducing equilibrium populations similar to EMD and experimental smFRET signatures.

Coarse-grained structure-based models (CG-SBMs; or G[o] models) are simplified potential energy functions of biomolecules or biomolecular complexes that encode their structure. Molecular dynamics simulations of such SBMs have been successfully used to study long time-scale dynamics such as protein and RNA folding, and large conformational transitions of biomolecular complexes. SBMs have several advantages: (1) Their MD simulations are computationally inexpensive, making extensive sampling easily accessible to many researchers. (2) They are easy to modify and can be adapted for the specific biomolecular problem that needs to be investigated. However, the force-fields of SBMs are not usually included in commonly used biomolecular simulation packages resulting in a barrier to their use. Here, we present SuBMIT (Structure Based Models Input Toolkit; https://github.com/sglabncbs/submit), a toolkit for generating coarse-grained SBM input files for performing MD simulations with GROMACS and OpenMM/OpenSMOG. Simulations whose input files can be generated using the different flavors of CG-SBMs present in SuBMIT include the folding and conformational ensembles of proteins with intrinsically disordered regions, 3D-domain-swapping in proteins and the dynamics of RNA-protein assemblies (e.g., simple RNA viruses).

Like other proteins, monoclonal antibodies - important biodrugs- are subject to post translational modifications, especially the N-glycosylations. However, the effect of the N-glycosylations remains poorly studied and atomistic details about their influence are rarely available. Moreover, the few existing studies focus on the prevalent immunoglobulin G1. To go further in the understanding of the impact of glycosylations, we have carried out a comparative exploration of the effect of N-glycosylations on two different classes of antibodies, namely Mab231, an IgG2 and the pembrolizumab, an IgG4. The two antibodies differ by their sequences, their length, their 3D structure but also by the location and composition of the glycans. In the present work, detailed and important information were gained through molecular dynamics simulations where both monoclonal antibodies were studied without and with the presence of their glycans. The results of 1.5 {micro}s of sampling for each system show that glycosylation does not drastically alter the overall conformational landscape of either antibody, whatever the metrics considered. However, it measurably modulates local flexibility, inter-domain correlated motions, and the relative orientation of the Fab arms with respect to the Fc domain, with statistically significant shifts in key geometric descriptors. Importantly, contact analysis reveals that glycan interactions extend beyond the Fc region to reach Fab residues. The allosteric network calculations demonstrate that the influence of Fc-bound glycans propagates even until the Fab framework regions in both mAbs, which could impact the antigen binding. The nature and magnitude of these effects are subclass-dependent, reflecting differences in glycan composition, hinge architecture, and three-dimensional organization Our findings challenge the prevailing view that Fc glycosylation uniformly promotes CH2 domain opening. More importantly, it underscores the necessity of considering full-length structures and IgG subclass diversity in glyco-engineering strategies.

Deep learning approaches have revolutionized protein structure prediction. These tools are trained using experimental data and recapitulate reported conformations, but there is great interest in predicting conformations that may be functionally relevant although experimentally underrepresented. Since many modern structure prediction tools use generative artificial intelligence diffusion models, we reframe the search for alternative molecular conformations as that of sampling from a diffusion distribution conditioned using any arbitrary Bayesian likelihood. We implement a twisted diffusion sampler in Boltz-2 to sample this conditioned distribution and demonstrate the utility of this approach, which does not require any additional training of the neural network, by implementing a diffusion analog of steered molecular dynamics simulations applied to mechanical systems. We can reproduce predicted stretched states of fragments of DNA, the muscle protein titin, and the inner-ear protocadherin-15 protein, as well as open states of the MscL ion channel consistent with experimental results. We expect that steered structure predictions will help sample underrepresented and non-equilibrium conformations for many macromolecular systems.

With the advance of hardware and software for molecular dynamics simulation it has become routine to obtain trajectories that are tens of microseconds in duration for all kinds of protein machinery. This shifts the burden of work onto analysis of the simulation data and opens opportunities for more rigorous and reproducible observations on mechanism. Toward this end we developed an investigator-blind analysis pipeline which operates on featurized simulation data, performs unsupervised clustering, and then identifies which input features are most discriminatory of cluster identity. Application of this pipeline to a large set of G-protein coupled receptor simulation data shows that it identifies several well-known microswitches. Inspection of these structural elements reveals changes in conformation that are known to accompany functional transitions of the receptor. In addition to these known structural elements the analysis also identifies two possibly new structural motifs: the kink in transmembrane helix 2, and a coupled "piston-like" motion of TM2 and TM3.

We present DIVIsive N-ary Ensembles (DIVINE), a deterministic, top-down clustering framework designed for molecular dynamics (MD) trajectories. DIVINE constructs a complete clustering hierarchy by recursively splitting clusters based on n-ary similarity principles, avoiding the need for O(N2) pairwise distance matrices. It supports multiple cluster selection criteria, including a weighted variance metric, and deterministic anchor initialization strategies such as NANI (N-ary Natural Initiation), ensuring reproducible and structurally meaningful partitions. Testing DIVINE up to a 305 s folding trajectory of the villin headpiece (HP35) revealed that it matched or exceeded the clustering quality of bisecting k-means while reducing runtime and eliminating stochastic variability. Its single-pass design enables efficient exploration of clustering resolutions without repeated executions. By combining scalability, interpretability, and determinism, DIVINE offers a robust and practical alternative to conventional MD clustering methods. DIVINE is publicly available as part of the MDANCE package: https://github.com/mqcomplab/MDANCE.

We developed AlphaUnfold, an automated pipeline that couples AF3 predictions with short-time (5 ns) high-pressure Molecular Dynamics (MD) using NAMD3. By subjecting models to baric stress, AlphaUnfold acts as a dynamic "stress-test" to identify structural fragility and potential unfolding. Testing a diverse set of proteins revealed a significant inverse correlation between average pLDDT and Root Mean Square Deviation (RMSD) after MD, indicating that lower confidence translates to rapid structural drift. Furthermore, domains with low local pLDDT consistently exhibited high Root Mean Square Fluctuation (RMSF), a behavior also observed in 200 ns simulations under standard pressure, pinpointing specific metastable areas. AlphaUnfold thus provides a viable, computationally efficient framework for assessing the biophysical robustness of AI-generated models, offering an "experimental-like" validation that ensures more reliable downstream applications in structural biology. MotivationAlphaFold3 (AF3) provides high-accuracy protein models characterized by the Predicted Local Distance Difference Test (pLDDT). However, these static predictions may harbor "not well-forged" regions lacking thermodynamic resilience. There is a critical need for rapid computational protocols to validate structural integrity beyond static confidence scores. AvailabilityGitHub: https://github.com/pegados/pipeline_AlphaUnfold Supplementary informationSupplementary data are available at http://biodados.icb.ufmg.br/alphaunfold Contacte-mail fabio, silva-jrp.miguel@ufmg.br

Molecular dynamics (MD) simulations are a powerful tool for investigating biomolecular dynamics underlying biological functions. However, the accessible spatiotemporal scales of conventional all-atom simulations remain limited by high computational costs. Coarse-graining reduces these costs by decreasing the number of interaction sites and enabling longer timesteps. In extreme cases, proteins are represented as single spherical particles; while such approximations facilitate cellular-scale simulations, they often sacrifice essential structural information, such as molecular shape and interaction anisotropy. Here, we present CGRig, a rigid-body protein model with residue-level interaction sites designed for long-time, large-scale simulations. In CGRig, each protein is treated as a single rigid-body embedding residue-level interaction sites. Its translational and rotational motions are described by the overdamped Langevin equation incorporating a shape-dependent friction matrix. Intermolecular interactions are calculated using G[o]-like native contact potentials, Debye-Huckel electrostatics, and volume exclusion. We validated that CGRig accurately reproduces the translational and rotational diffusion coefficients expected from the friction matrix for an isolated protein. For dimeric systems, the model successfully maintained native complex structures. Furthermore, two initially separated proteins converged into the correct complex with an association rate consistent with all-atom simulations. Notably, CGRig achieved a simulation performance exceeding 17 s/day for a 1,024-molecule system. These results demonstrate that CGRig provides an efficient framework for simulating protein assembly while retaining residue-level interaction specificity, making it a valuable tool for investigating large-scale biomolecular self-assembly.

The integration of coarse-grained (CG) approaches into docking workflows offers a powerful strategy for modelling large biomolecular assemblies with reduced computational costs. We present here the implementation of the MARTINI2 coarse-grained force field into the HADDOCK3 integrative modelling platform. This development enables the use of the CG representations and parameters within HADDOCK3 for efficient sampling and scoring of large protein-protein complexes. The implementation takes advantage of the modular and flexible architecture of HADDOCK3, allowing a seamless combination of MARTINI2 representation with the various modules. Conversion from and to all-atom models is integrated into the coarse-grained modelling workflow. The performance of the protocol is first assessed on protein-protein and protein-DNA benchmarks and then illustrated on a few representative large-scale systems, demonstrating a significant reduction in computational costs while maintaining biologically relevant accuracy.

Phase separated condensates are recognized as a ubiquitous mechanism of spatial organization in cell biology. Biophysical modeling of condensates provides critical insights into the dynamics and functions of these subcellular structures that are difficult to extract via experiments. Here we present an efficient computational pipeline, CASPULE (Condensate Analysis of Sticker Spacer Polymers Using the LAMMPS Engine), to simulate and analyze the biological condensates made of sticker-spacer polymers. CASPULE implements a unique force field that combines traditional Langevin dynamics with a "detailed balance proof" protocol for single-valent bond formation between stickers. This framework allows us to study the non-trivial biophysics that emerge out of the single-valent sticker interactions coupled with the effect of separation in energetic contribution by stickers and spacers. We provide detailed documentation on how to setup the simulation environment, perform simulations and analyze the results. Through case studies, we highlight the utility and efficacy of our pipeline. Importantly, we provide statistical parameters to characterize the cluster size distribution often observed in biological systems. We envision this tool to be broadly useful in decoding the interplay of kinetics and thermodynamics underlying the formation and function of biological condensates.

Characterizing conformational transitions between distinct structural states is essential for understanding protein function but remains challenging due to the timescale limitations of atomistic molecular dynamics. While coarse-grained models like Martini accelerate sampling, classical elastic-network or G[o]-like restraints often trap proteins in a single energy basin, precluding the study of transition pathways between distinct functional states. Here, we present CTGoMartini, a comprehensive Python package designed to simulate protein conformational transitions using G[o]-Martini models in explicit membranes. CTGoMartini addresses key methodological limitations of existing approaches by redefining native contacts as a dedicated interaction type, thereby eliminating spurious protein aggregation artifacts in multi-copy simulations. The package implements both switching and multiple-basin approaches (Exponential and Hamiltonian mixing) to sample transitions between experimentally defined states. Furthermore, it integrates Hamiltonian replica exchange molecular dynamics (HREMD) with PyMBAR analysis, enabling efficient optimization of mixing parameters that govern barrier heights and relative state stabilities. We demonstrate the power of CTGoMartini through two biologically significant membrane protein systems: (1) capturing the inward-open to outward-open transition of the lipid transporter SPNS2, revealing the molecular mechanism of S1P translocation; and (2) elucidating how membrane surface tension and anionic lipids (POPA, PIP2) modulate the conformational equilibrium of the mechanosensitive ion channel TREK1. By streamlining model construction, simulation, and analysis, CTGoMartini offers an easy-to-use platform that connects static structural snapshots with their underlying dynamic functional mechanisms. TOC Graphic O_FIG O_LINKSMALLFIG WIDTH=200 HEIGHT=118 SRC="FIGDIR/small/721921v1_ufig1.gif" ALT="Figure 1"> View larger version (26K): org.highwire.dtl.DTLVardef@75eb26org.highwire.dtl.DTLVardef@1a12accorg.highwire.dtl.DTLVardef@e927org.highwire.dtl.DTLVardef@1cb0dcd_HPS_FORMAT_FIGEXP M_FIG C_FIG

Reconstructing all-atom (AA) structures from highly coarse-grained (HCG) models remains a significant challenge in multiscale molecular dynamics (MD) simulations, particularly for mesoscale biomolecular assemblies that are beyond the reach of conventional MD methods. Building upon ProNet Backmapping, a neural-network-based thermodynamically consistent approach, we introduce a progressive backmapping framework that reconstructs AA models in a stepwise manner across neighboring resolutions, for example, from a 3-residue-per-site HCG model to a 1-residue-per-site model, then to an AA model. This progressive backmapping method achieves high accuracy across a wide range of proteins and effectively reconstructs flexible linkers in multidomain architectures. Moreover, it supports hierarchical reconstruction of complex protein assemblies, including multiple virus-like particles spanning tens of nanometers and containing hundreds of subunits. Using this framework, we demonstrate--for the first time--the ability to hierarchically backmap entire viral assemblies from HCG to full AA resolution, covering at least three different resolutions. Overall, our method provides a scalable framework for incorporating atomistic detail into mesoscale simulations of complex systems across many applications in chemistry and biology. Table of contents figure O_FIG O_LINKSMALLFIG WIDTH=200 HEIGHT=70 SRC="FIGDIR/small/709104v1_ufig1.gif" ALT="Figure 1"> View larger version (45K): org.highwire.dtl.DTLVardef@4af423org.highwire.dtl.DTLVardef@e2669borg.highwire.dtl.DTLVardef@1be80eforg.highwire.dtl.DTLVardef@2e679_HPS_FORMAT_FIGEXP M_FIG C_FIG

AO_SCPLOWBSTRACTC_SCPLOWWe present a general approach to find amino acid sequences corresponding to the most compact enzyme likely to retain the structure of a given catalytic site. Our approach is based on using Monte Carlo (MC) simulations to sample an energy landscape where minima correspond, by construction, to sequences with the aforementioned properties. Building on previous work (Wu et al., 2025) and with the BAGEL package (Lala et al., 2025), we implement a route to achieve this goal using only the information extracted from a protein language model (PLM), without structural information. After generating a set of candidate sequences with this PLM-guided BAGEL optimization, we further filter potential candidates for downstream experimental validation using a two-stage protocol. First, deep-learning-based structure prediction models (ESMFold, Chai-1, Boltz-2) are used to identify a structural consensus among designs with highly conserved active-site geometries, yielding many candidates with active-site RMSD below a few angstroms relative to the wild-type and pLDDT scores above 80. Second, molecular dynamics simulations are performed on a filtered subset of sequences (based on active-site RMSD and SolubleMPNN log-likelihoods) to evaluate active-site stability when including thermal fluctuations. For the most promising enzymes, these yield RMSF values in the active site below 1.0 [A] and an active-site RMSD drift between 0.5 and 1.5 [A], making these mini-variants comparable to the wild type, though outcomes vary across enzymes. Given the protocols generality, we believe these results represent a step forward in AI-guided enzyme design. To facilitate rapid experimental validation by the broader community, we open-source all sequences generated by our computational pipeline. These include designs for four representative enzymes of this study: PETase, subtilisin Carlsberg (serine protease), Taq DNA polymerase, and VioA.

Accurate ion force fields are essential for molecular dynamics simulations of biomolecular systems, particularly in combination with modern water models such as OPC. While OPC water improves the description of bulk water and biomolecules, the transferability of existing ion force fields to this model remains an open question. Here, we systematically assess the transferability of monovalent and divalent ion force field parameters (Li+, Na+, K+, Cs+, Mg2+,Ca2+, Sr2+, Ba2+, Cl- and Br-) to OPC water by comparing single-ion and ion-pairing properties with experimental data. Our analysis reveals that no single literature parameter set provides accurate results for all ions when directly transferred to OPC water. We hence introduce the MS/G-LB(OPC) force field, which combines Mamatkulov-Schwierz-Grotz cation parameters with Loche-Bonthuis anion parameters. MS/G-LB(OPC) reproduces hydration free energies, first-shell structural properties and activity derivatives at low salt concentrations. Our results demonstrate that transferring ion parameters to OPC can lead to significant and ion-specific deviations from experimental data, making careful validation essential. At the same time, the systematic transfer and combination of ion parameters from existing force fields can provide a practical and computationally efficient alternative to full reparameterization. MS/G-LB(OPC) is available at https://git.rz.uni-augsburg.de/cbio-gitpub/opc-ion-force-fields.

Hyperpolarization-activated cyclic nucleotide-gated (HCN) channels are are a family of voltage-gated, cyclic-nucleotide modulated Na+/K+ channels that regulate spontaneous rhythmic electrical activity in both the heart and the brain. Understanding differences in the responsiveness to cyclic adenosine monophosphate (cAMP) modulation between HCN isoforms would offer insight into the specific binding interactions that drive channel activation. Using all-atom molecular dynamics (MD) simulations and the free-energy perturbation (FEP) approach, we determined the absolute binding free energy of cAMP to the the cyclicnucleotide-binding domain (CNBD) of HCN isoforms 1-4. By studying the free-energy of ligand binding to the various isoforms of HCN, our study advances the understanding of HCN channel activation and modulation mechanisms. Overall, our work offers insight into explaining differences in channel sensitivity across the isoforms of HCN.

Structure-based drug design has become increasingly important in the pharmaceutical industry for accelerating the discovery of effective drug candidates. In particular, ligand binding free energy serves as a critical metric for predicting drug efficacy during the key stages of hit discovery and lead optimization. Continuous progresses have been made in the prediction of ligand binding free energies, but direct comparisons of different methods using the same force field remain challenging due to their unique implementations into different simulation engines. In this study, we present a direct comparison of four popular methodologies (Amber-TI, GROMACS-NETI, OpenMM-FEP, and BLaDE-MSLD) for calculating relative binding free energies ({Delta}{Delta}Gbind) with the same Amber protein and ligand force fields using MolCube Alchemical Free Energy Simulator (MolCube-AFES), which provides an input generation workflow to support {Delta}{Delta}Gbind calculations of all four methods. We used 80 alchemical transformations (among the JACS benchmark set by Wang et al.) and two additional applications to compare the predicted {Delta}{Delta}Gbind from the four methods against experimental measurements. All four methods reproduced experimentally observed trends with most transformations within {+/-}2 kcal/mol from experiments and show broadly comparable accuracy with no statistically significant performance differences across the benchmark dataset. These results demonstrate that MolCube-AFES enables controlled, cross platform benchmarking and show that all four different alchemical free energy methods deliver statistically equivalent accuracy, with method selection guided by workflow requirements such as throughput, portability, and perturbation network design rather than expected differences in performances.